## Understanding Public Attitudes toward Restrictive Voting Laws
## Katie Clayton  
## Step 2: Analyze data
## Last updated: June 8, 2023

# Initial settings --------------------------------------------------------

rm(list = ls())
library(tidyverse)
library(estimatr)
library(ggthemes)
library(magrittr)
library(margins)
library(multcomp)

# Read data ---------------------------------------------------------------

df <- read.csv("output/cleaned_dataset.csv") %>%
  mutate(group = factor(group, levels = c("control","harm","boost"))) 

# Run models  -------------------------------------------------------------

mod1 <- estimatr::lm_robust(support_law~group*rep.dum+b1+b2+b3+b5+b6+b7,df)
me1.dem <- summary(margins(mod1, data = subset(df, rep.dum == 0)))
me1.rep <- summary(margins(mod1, data = subset(df, rep.dum == 1)))

mod2 <- estimatr::lm_robust(election.trust~group*rep.dum+b1+b2+b3+b5+b6+b7,df)
me2.dem <- summary(margins(mod2, data = subset(df, rep.dum == 0)))
me2.rep <- summary(margins(mod2, data = subset(df, rep.dum == 1)))

mod3 <- estimatr::lm_robust(emo.anger~group*rep.dum+b1+b2+b3+b5+b6+b7,df)
me3.dem <- summary(margins(mod3, data = subset(df, rep.dum == 0)))
me3.rep <- summary(margins(mod3, data = subset(df, rep.dum == 1)))

mod4 <- estimatr::lm_robust(lottery_3~group*rep.dum+b1+b2+b3+b5+b6+b7,df)
me4.dem <- summary(margins(mod4, data = subset(df, rep.dum == 0)))
me4.rep <- summary(margins(mod4, data = subset(df, rep.dum == 1)))


# Visualize main results --------------------------------------------------

covar <- c("Democrat","Democrat","Republican","Republican","Democrat","Democrat","Republican","Republican","Democrat","Democrat","Republican","Republican")

term <- rep("condition", 12)

contrast <- rep(c("Harm minorities","Boost minority turnout"), 6)

null.value <- rep(0, 12)

estimate <- c(me1.dem$AME[[8]],me1.dem$AME[[7]],me1.rep$AME[[8]],me1.rep$AME[[7]],
              me2.dem$AME[[8]],me2.dem$AME[[7]],me2.rep$AME[[8]],me2.rep$AME[[7]],
              me3.dem$AME[[8]],me3.dem$AME[[7]],me3.rep$AME[[8]],me3.rep$AME[[7]])

std.error <- c(me1.dem$SE[[8]],me1.dem$SE[[7]],me1.rep$SE[[8]],me1.rep$SE[[7]],
               me2.dem$SE[[8]],me2.dem$SE[[7]],me2.rep$SE[[8]],me2.rep$SE[[7]],
               me3.dem$SE[[8]],me3.dem$SE[[7]],me3.rep$SE[[8]],me3.rep$SE[[7]])

conf.low <- c(me1.dem$lower[[8]],me1.dem$lower[[7]],me1.rep$lower[[8]],me1.rep$lower[[7]],
              me2.dem$lower[[8]],me2.dem$lower[[7]],me2.rep$lower[[8]],me2.rep$lower[[7]],
              me3.dem$lower[[8]],me3.dem$lower[[7]],me3.rep$lower[[8]],me3.rep$lower[[7]])

conf.high <- c(me1.dem$upper[[8]],me1.dem$upper[[7]],me1.rep$upper[[8]],me1.rep$upper[[7]],
               me2.dem$upper[[8]],me2.dem$upper[[7]],me2.rep$upper[[8]],me2.rep$upper[[7]],
               me3.dem$upper[[8]],me3.dem$upper[[7]],me3.rep$upper[[8]],me3.rep$upper[[7]])

statistic <- c(me1.dem$z[[8]],me1.dem$z[[7]],me1.rep$z[[8]],me1.rep$z[[7]],
               me2.dem$z[[8]],me2.dem$z[[7]],me2.rep$z[[8]],me2.rep$z[[7]],
               me3.dem$z[[8]],me3.dem$z[[7]],me3.rep$z[[8]],me3.rep$z[[7]])

p.value <- c(me1.dem$p[[8]],me1.dem$p[[7]],me1.rep$p[[8]],me1.rep$p[[7]],
             me2.dem$p[[8]],me2.dem$p[[7]],me2.rep$p[[8]],me2.rep$p[[7]],
             me3.dem$p[[8]],me3.dem$p[[7]],me3.rep$p[[8]],me3.rep$p[[7]])

dv_lab <- c(rep("Support law", 4), rep("Election faith", 4), rep("Anger", 4))

subg <- c("Democrats","Democrats","Republicans","Republicans","Democrats","Democrats","Republicans","Republicans","Democrats","Democrats","Republicans","Republicans")

mod <- c(rep("mod1", 4), rep("mod2", 4), rep("mod3", 4))

df.viz <- cbind(covar, term, contrast, null.value, estimate, std.error, conf.low, conf.high, statistic, p.value, dv_lab, subg, mod) %>%
  as.data.frame() %>% 
  mutate(estimate = as.numeric(estimate)) %>% 
  mutate(lab = ifelse(abs(estimate) < .01, round(estimate, digits = 3), round(estimate, digits = 2))) %>% 
  mutate(lab = as.character(lab)) %>% 
  mutate(lab = str_replace(lab, fixed("0."), ".")) %>% 
  mutate(lab = str_c(lab, 
                     c("***", "**", "*", "") %>% 
                       magrittr::extract(
                         p.value %>% 
                           findInterval(
                             c(-Inf, .001, .01, .05, Inf)
                           )
                       ) )) %>% 
  mutate(dv_lab = factor(dv_lab, c("Anger","Election faith","Support law"))) %>% 
  mutate(std.error = as.numeric(std.error),
         conf.low = as.numeric(conf.low),
         conf.high = as.numeric(conf.high))

viz <- df.viz %>% 
  mutate(contrast = contrast %>% 
           factor(
             c("Harm minorities", "Boost minority turnout"))) %>% 
  ggplot() +
  geom_vline(
    aes(xintercept = 0),
    linetype = "dashed") + 
  geom_linerange(
    aes(xmin = conf.low, 
        xmax = conf.high, 
        y = dv_lab,
        color = covar,
        group = covar),
    data = df.viz,
    position = position_dodge2(width = .4), 
    show.legend = F
  ) +
  geom_point(
    aes(estimate, dv_lab,
        color = covar,
        group = covar),
    data = df.viz,
    position = position_dodge2(width = .4),
    shape = 21,
    size = 8.75,
    fill = "white"
  ) +
  geom_point(
    aes(estimate, dv_lab,
        group = covar),
    data = df.viz,
    position = position_dodge2(width = .4),
    shape = 21,
    size = 7.75,
    color = "white",
    fill = "white", 
    show.legend = F
  ) +
  geom_text(
    aes(
      estimate, dv_lab, label = lab,
      color = covar,
      group = covar
    ),
    data = df.viz,
    position = position_dodge2(width = .4),
    size = 2,
    show.legend = F
  ) +
  scale_y_discrete(
    breaks = c(
      "Anger",
      "Election faith",
      "Support law"
    ),
    labels = c(
      "Anger toward law",
      "Faith in future elections",
      "Support strict voting law"
    ) %>%
      str_wrap(width = 20)
  ) +
  scale_color_manual(
    values = c("#e91b23", "#2fabe1") %>% rev,
    guide = guide_legend(
      override.aes = list(shape = 16, size = 3),
      reverse = T
    )
  ) +
  facet_grid(
    contrast~., 
    scales = "free_y",
    space = "free_y",
    labeller = label_wrap_gen(width = 13)
  ) +
  labs(
    title = NULL,
    x = "Difference relative to control",
    y = "",
    color = ""
  ) +
  theme(
    legend.position = "bottom",
    panel.background  = 
      element_rect(color = "grey95",
                   fill = "grey95"),
    panel.grid = element_blank(),
    strip.background = element_rect(color = "white",
                                    fill = "white"),
    legend.background = element_rect(color = "white",
                                     fill = "white"),
    panel.grid.minor = element_blank(),
    plot.title = element_text(hjust = 0.5)
  ) 

viz

ggsave("figures/main-fig.png", viz, width = 10, height = 6)


# Visualize lottery results -----------------------------------------------

lottery.diffmeans <- df %>% 
  group_by(group, rep.dum) %>%   
  summarise(lottery_1m = mean(lottery_1, na.rm = TRUE), 
            lottery_2m = mean(lottery_2, na.rm = TRUE), 
            lottery_3m = mean(lottery_3, na.rm = TRUE), 
            n = n(),
            lottery_1se = sd(lottery_1, na.rm = TRUE)/sqrt(n()),
            lottery_2se = sd(lottery_2, na.rm = TRUE)/sqrt(n()),
            lottery_3se = sd(lottery_3, na.rm = TRUE)/sqrt(n())) %>% 
  mutate(lottery_1lower = lottery_1m - qt(1 - ((1 - 0.95) / 2), n - 1) * lottery_1se) %>% 
  mutate(lottery_1upper = lottery_1m + qt(1 - ((1 - 0.95) / 2), n - 1) * lottery_1se) %>% 
  mutate(lottery_2lower = lottery_2m - qt(1 - ((1 - 0.95) / 2), n - 1) * lottery_2se) %>% 
  mutate(lottery_2upper = lottery_2m + qt(1 - ((1 - 0.95) / 2), n - 1) * lottery_2se) %>% 
  mutate(lottery_3lower = lottery_3m - qt(1 - ((1 - 0.95) / 2), n - 1) * lottery_3se) %>% 
  mutate(lottery_3upper = lottery_3m + qt(1 - ((1 - 0.95) / 2), n - 1) * lottery_3se) %>% 
  mutate(rep.dum = ifelse(rep.dum == 1, "Republican", "Democrat")) %>% 
  mutate(rep.dum = factor(rep.dum, levels = c("Republican", "Democrat"), labels = c("Republican (n = 1328)", "Democrat (n = 1773)")))


viz2 <- ggplot(lottery.diffmeans, aes(x = group, y = lottery_3m)) + 
  geom_col(aes(fill = rep.dum), alpha = 0.7) +
  scale_fill_manual(values = c("indianred1","skyblue1")) +
  scale_color_manual(values = c("darkred","navy")) +
  facet_grid(~rep.dum) +
  geom_linerange(aes(ymin = lottery_3lower, ymax = lottery_3upper, color = rep.dum)) +
  scale_x_discrete(labels = c("control" = "Control",
                              "harm" = "Harm minorities",
                              "boost" = "Boost minority turnout")) +
  labs(title = NULL, 
       x = NULL,
       y = "Amount in dollars ($100 maximum)") +
  theme_few() + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  ylim(0, 30)

viz2

viz3 <- ggplot(lottery.diffmeans, aes(x = group, y = lottery_2m)) + 
  geom_col(aes(fill = rep.dum), alpha = 0.7) +
  scale_fill_manual(values = c("indianred1","skyblue1")) +
  scale_color_manual(values = c("darkred","navy")) +
  facet_grid(~rep.dum) +
  geom_linerange(aes(ymin = lottery_2lower, ymax = lottery_2upper, color = rep.dum)) +
  scale_x_discrete(labels = c("control" = "Control",
                              "harm" = "Harm minorities",
                              "boost" = "Boost minority turnout")) +
  labs(title = NULL, 
       x = NULL,
       y = "Amount in dollars ($100 maximum)") +
  theme_few() + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  ylim(0, 30)

viz3

ggsave("figures/lottery.png", viz2, width = 10, height = 6)
ggsave("figures/lottery_v2.png", viz3, width = 10, height = 6)



# Create Table 1 ----------------------------------------------------------

texreg::texreg(mod1, digits = 3)
texreg::texreg(mod2, digits = 3)
texreg::texreg(mod3, digits = 3)

me1.dem
me2.dem
me3.dem

me1.rep
me2.rep
me3.rep

lincom1.party1 <- summary(glht(mod1, linfct = c("groupharm:rep.dum + groupharm = 0")))
lincom1.party2 <- summary(glht(mod1, linfct = c("groupboost:rep.dum + groupboost = 0")))
list1.party <- c(coef(lincom1.party1), sqrt(vcov(lincom1.party1)), coef(lincom1.party2), sqrt(vcov(lincom1.party2)))
list1.party

lincom2.party1 <- summary(glht(mod2, linfct = c("groupharm:rep.dum + groupharm = 0")))
lincom2.party2 <- summary(glht(mod2, linfct = c("groupboost:rep.dum + groupboost = 0")))
list2.party <- c(coef(lincom2.party1), sqrt(vcov(lincom2.party1)), coef(lincom2.party2), sqrt(vcov(lincom2.party2)))
list2.party

lincom3.party1 <- summary(glht(mod3, linfct = c("groupharm:rep.dum + groupharm = 0")))
lincom3.party2 <- summary(glht(mod3, linfct = c("groupboost:rep.dum + groupboost = 0")))
list3.party <- c(coef(lincom3.party1), sqrt(vcov(lincom3.party1)), coef(lincom3.party2), sqrt(vcov(lincom3.party2)))
list3.party

